library(foreign)
library(HonestDiD)
library(lfe)
setwd("C:/Users/dpowell/Documents/mscontin/")
mydata <- read.dta("forr_wrongstates.dta")


for (j in 1983:1994) {
  nam <- paste("cc", j,sep="")
  assign(nam, as.integer(mydata$year==j) * mydata$nontripl)
}

for (j in 1996:2017) {
  nam <- paste("cc", j,sep="")
  assign(nam, as.integer(mydata$year==j) * mydata$nontripl)
}
mort = lfe::felm(overdose_rate ~ cc1983+cc1984+cc1985 +cc1986 + cc1987 + cc1988 +cc1989 + cc1990 + cc1991 +cc1992 + cc1993 +cc1994 + cc1996  +cc1997 +cc1998 + cc1999 +cc2000 +cc2001 +cc2002 + cc2003 + cc2004 + cc2005 + cc2006 + cc2007 + cc2008 +cc2009 + cc2010 + cc2011 + cc2012 + cc2013 + cc2014 + cc2015 + cc2016 + cc2017 | factor(year) + factor(stfips) | 0 | stfips  , data=mydata, weights=mydata$totpop)

coefIndex = which(grepl(x = dimnames(mort$coefficients)[[1]],
                        pattern = "cc"))
betahat = mort$beta[coefIndex, ]
sigma = mort$clustervcv[coefIndex, coefIndex]


timeVec = c(seq(from = -13, to = -2, by = 1), seq(from = 0, to = 21, by = 1))
referencePeriod = -1
postPeriodIndices = which(timeVec > -1)
prePeriodIndices = which(timeVec < -1)
stdErrors = summary(mort)$coefficients[coefIndex,2]

LWdata_EventStudy = list(
  betahat = betahat,
  sigma = sigma,
  timeVec = timeVec,
  referencePeriod = referencePeriod,
  prePeriodIndices = prePeriodIndices,
  postPeriodIndices = postPeriodIndices,
  stdErrors = stdErrors
)


l_vec =  (basisVector(6, 22) + basisVector(7, 22) + basisVector(8, 22) + basisVector(9, 22) + basisVector(10, 22) + basisVector(11, 22) + basisVector(12, 22) + basisVector(13, 22) + basisVector(14, 22) + basisVector(15, 22))/10
l_vec =  (basisVector(6, 22) + basisVector(7, 22) + basisVector(8, 22) + basisVector(9, 22) + basisVector(10, 22) + basisVector(11, 22) + basisVector(12, 22) + basisVector(13, 22) + basisVector(14, 22) + basisVector(15, 22))/10
l_vec = (basisVector(16, 22) + basisVector(17, 22) + basisVector(18, 22) + basisVector(19, 22) + basisVector(20, 22) + basisVector(21, 22) + basisVector(22, 22))/7
l_vec =  (basisVector(6, 22) + basisVector(7, 22) + basisVector(8, 22) + basisVector(9, 22) + basisVector(10, 22) + basisVector(11, 22) + basisVector(12, 22) + basisVector(13, 22) + basisVector(14, 22) + basisVector(15, 22))/10
l_vec = (basisVector(1, 22) + basisVector(2, 22) + basisVector(3, 22) + basisVector(4, 22) + basisVector(5, 22))/5

l_vec
computeConditionalCS_DeltaRMI(betahat=betahat, sigma=sigma, numPrePeriods=12, numPostPeriods=22,
                              l_vec = l_vec, Mbar = 0,
                              alpha = 0.05, hybrid_flag = "LF", hybrid_kappa = 0.05/10,
                              returnLength = F, postPeriodMomentsOnly = T,
                              gridPoints=10^3, grid.ub = NA, grid.lb = NA)

DeltaSD_RobustResults = createSensitivityResults(betahat = betahat,
                                                 sigma = sigma,
                                                 numPrePeriods = 12,
                                                 numPostPeriods = 22,
                                                 l_vec = l_vec,
                                                 Mvec = seq(from = 0, to = 0.04, by = 0.005))
head(DeltaSD_RobustResults)

OriginalResults = constructOriginalCS(betahat = betahat,
                                      sigma = sigma,
                                      numPrePeriods = 12,
                                      numPostPeriods = 22,
                                      l_vec = l_vec )
# Construct sensitivity plot
DeltaSD_SensitivityPlot = createSensitivityPlot(robustResults = DeltaSD_RobustResults,
                                                originalResults = OriginalResults)
DeltaSD_SensitivityPlot
DeltaSD_RobustResults
OriginalResults

